其他
pysam 读取 bam 文件准备 Ribo-seq 质控数据
点击上方关注老俊俊
1引言
前面关于 Ribo-seq 的质控我们是直接在 R 里面读取 bam 文件进行操作,我们也可以使用 pysam 处理好数据再导入 R 里面进行分析也可。
2读取 bam 文件
数据还是比对到转录组上:
$ zless 4.rmtRNA-data/ribo-ESC_rmtRNA.fq.gz |\
bowtie ../../index/trans-index/longest_cds_trans \
-q -p 10 -m 1 -v 2 \
--best --strata - -S |\
samtools sort -@ 10 \
-o 7.map2trans-data/ribo-ESC.sorted.bam
安装 pysam:
linux 系统里面才能装
,windows 用户这步就不用考虑了。
!pip install pysam
import pysam
import pandas as pd
读取文件:
samfile = pysam.AlignmentFile("ribo-ESC.sorted.bam", "rb")
info = []
for line in samfile:
rname = str(line.reference_name)
if rname == 'None':
continue
# split rname
fileds = rname.split('|')
cds_start = int(fileds[3])
cds_end = int(fileds[4])
# shift aln pos
pos = int(line.pos + 1)
# assign frame
mod = abs(pos - cds_start) % 3
if mod == 0:
frame = 'frame 0'
elif mod == 1:
frame = 'frame 1'
elif mod == 2:
frame = 'frame 2'
else:
pass
# calculate relative to cds_start/stop distance
dist2start = pos - cds_start
dist2stop = pos - cds_end
# read length
readlen = line.qlen
# save in list
info.append([rname,readlen,frame,dist2start,dist2stop])
这里计算了每个
read
的 frame 类型,长度,与 start codon 和 stop codon 的相对距离。
注意 pysam 读进来的
比对的位置会减 1
,所以位置得再加上 1 才是原来 bam 文件的位置。
转为数据框:
df = pd.DataFrame(info,columns=['rname','readlen','frame','dist2start','dist2stop'])
# save
# df.to_csv('df_res.csv',index=False,header=True)
df
导出来在 R 里面就可以分析了。
3读取 sam 文件
如果想在 windows 里操作,你可以直接将比对结果保存为 sam 文件:
$ zless 4.rmtRNA-data/ribo-ESC_rmtRNA.fq.gz |\
bowtie ../../index/trans-index/longest_cds_trans \
-q -p 10 -m 1 -v 2 \
--best --strata - -S |\
samtools sort -@ 10 -O SAM \
-o 7.map2trans-data/ribo-ESC.sorted.sam
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from collections import Counter
处理 sam 文件:
info = []
with open('ribo-ESC.sorted.sam','r') as sam:
for line in sam:
if line.startswith('@'):
continue
fields = line.split()
if fields[2] == '*':
continue
# split line
alnpos = int(fields[3])
startpos = int(fields[2].split('|')[3])
stoppos = int(fields[2].split('|')[4])
# assign frame
mod = abs(alnpos - startpos) % 3
if mod == 0:
frame = 'frame 0'
elif mod == 1:
frame = 'frame 1'
elif mod == 2:
frame = 'frame 2'
else:
pass
# calculate relative to cds_start/stop distance
dist2start = alnpos - startpos
dist2stop = alnpos - stoppos
# read length
rlen = len(fields[9])
info.append([fields[2],rlen,frame,dist2start,dist2stop])
转为数据框:
df = pd.DataFrame(info,columns=['rname','rlen','frame','dist2start','dist2stop'])
# save
# df.to_csv('df_res.csv',index=False,header=True)
df
我们可以简单统计一下长度:
tmp_freq = Counter(df.rlen).most_common()
# sort ascending
tmp_freq.sort(key= lambda item:item[0])
# to data frame and save
tmp_freq_df = pd.DataFrame(tmp_freq,columns=['readlen','numbers'])
tmp_freq_df
可视化一下:
sns.set_theme(style="white", context="talk")
# Generate some sequential data
sns.barplot(x=tmp_freq_df.readlen, y=tmp_freq_df.numbers, palette="Paired")
欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群
哦,数据代码已上传至QQ群,欢迎加入下载。
群二维码:
老俊俊微信:
知识星球:
所以今天你学习了吗?
今天的分享就到这里了,敬请期待下一篇!
最后欢迎大家分享转发,您的点赞是对我的鼓励和肯定!
如果觉得对您帮助很大,赏杯快乐水喝喝吧!
往期回顾
◀跟着 Cell Reports 学做图-CLIP-seq 数据可视化
◀m6A enrichment on peak center
◀m6A metagene distribution 纠正坐标
◀...